# CLASS 1 - SCRIPT #1
# (C) Fernando Ascensão 2025


library(rgbif)
library(sf)
library(dplyr)
library(leaflet)
library(ggplot2)
library(tidyterra)
library(terra)
library(gridExtra)
library(eks) # kernel smoothing functions
library(gstat)



# install.packages("tidyterra", dep=T) [all packages]


rm(list=ls())


# load data ---------------------------------------------------------------

mywd <- "C:/Users/fjascensao/OneDrive - Universidade de Lisboa/Aulas/MIBC/Classes/Class_1_20260302"  # <- change to your root!

setwd(mywd)

# study area
Iberia <- st_read(paste0(mywd, "/Administrative/Iberia.shp"))
Iberia

Iberia4326 <- st_transform(Iberia, 4326)

plot(st_geometry(Iberia4326))



# retrieve data from GBIF -------------------------------------------------

# select species name
taxon_name_query <- "Pelodytes ibericus"

species_key_ <- rgbif::name_backbone(taxon_name_query)
species_key <- species_key_[[1]]

occ = rgbif::occ_data(taxonKey=species_key,
                      country = c("PT","ES"),
                      hasCoordinate = TRUE, 
                      limit = 10000)

#coordinateUncertaintyInMeters = "0,500",
# basisOfRecord = "HUMAN_OBSERVATION",
#year = '2000,2025')


str(occ)
nrow(occ$PT$data)
nrow(occ$ES$data)

occ_data1 = occ$PT$data
occ_data2 = occ$ES$data

occ_data <- bind_rows(occ_data1, occ_data2)

names(occ_data)
nrow(occ_data)

table(occ_data$basisOfRecord)
table(occ_data$country)



# select columns
(occ_myd <-  occ_data %>%
  select(order, familyKey, scientificName, decimalLongitude, decimalLatitude, 
         year, month, day, eventDate))


# plot by year
#
ggplot(occ_myd, aes(year)) +
  geom_bar() +
  theme_minimal() 



# convert to spatial
occ_sf <- st_as_sf(occ_myd,  
                   coords = c("decimalLongitude","decimalLatitude"), 
                   crs=4326)

names(occ_sf)


leaflet() %>% 
  setView(lng = -4, lat = 39.5, zoom = 6) %>%
  addTiles() %>%
  addCircles(lng = ~decimalLongitude, lat = ~decimalLatitude, 
             data = occ_myd)


# Repeat for "Salamandra salamandra"



# rasterize ---------------------------------------------------------------


occ_sf # <-- WGS84


mystudyarea3763 <- Iberia %>%
  filter(shapeISO  == "PRT") %>%
  st_union() %>%
  st_transform(., 3763)

plot(st_geometry(mystudyarea3763))

index <- occ_sf %>%
  st_transform(., 3763) %>%
  st_intersects(mystudyarea3763, sparse = F)

occ_sf3763 <- occ_sf %>%
  st_transform(., 3763)
  
occ_sf3763 <- occ_sf3763[index,]


(g1 <- ggplot(mystudyarea3763) +
    geom_sf(fill="grey99", col="grey") +
    geom_sf(data=occ_sf3763, col="blue") +
    theme(legend.position = "none"))


my_template1k <- rast(ext(occ_sf3763), res=1000)
my_template5k <- rast(ext(occ_sf3763), res=5000)
my_template10k <- rast(ext(occ_sf3763), res=10000)

myr1k <- rasterize(occ_sf3763, my_template1k)
myr5k <- rasterize(occ_sf3763, my_template5k)
myr10k <- rasterize(occ_sf3763, my_template10k)
myr10k_ <- rasterize(occ_sf3763, my_template10k, fun=sum)


(r1 <- ggplot(mystudyarea3763) +
  geom_sf(fill=NA, col="grey") +
  theme_void() +
  geom_spatraster(data = myr1k) +
  scale_fill_gradient(na.value = NA) +
  theme(legend.position = "none"))

(r2 <- ggplot(mystudyarea3763) +
  geom_sf(fill=NA, col="grey") +
  theme_void() +
  geom_spatraster(data = myr5k) +
  scale_fill_gradient(na.value = NA) +
  theme(legend.position = "none"))

(r3 <- ggplot(mystudyarea3763) +
  geom_sf(fill=NA, col="grey") +
  theme_void() +
  geom_spatraster(data = myr10k) +
  scale_fill_gradient(na.value = NA) +
  theme(legend.position = "none"))

(r4 <- ggplot(mystudyarea3763) +
  geom_sf(fill=NA, col="grey") +
  theme_void() +
  geom_spatraster(data = myr10k_) +  # log?
  scale_fill_viridis_c(na.value = NA) )

grid.arrange(r1, r2, 
             r3, r4)


# save raster
# names(myr10k_) <- "density10k"
# writeRaster(myr10k_, "myr10k_.tiff", overwrite=T)



# density -----------------------------------------------------------------

myr10k_df <- as.data.frame(myr10k_, xy = TRUE)
head(myr10k_df)

names(myr10k_df) <- c("x", "y", "z")

idw_r <- interpIDW(myr10k_, as.matrix(myr10k_df), radius=50000, power=1, smooth=1, maxPoints=5)
crs(idw_r) <- crs(vect(mystudyarea3763))

# Create a template raster covering the full study area extent
mystudyarea3763_r <- rast(vect(mystudyarea3763), resolution = 1000)

# Rasterize the study area polygon (1 inside, NA outside)
mystudyarea3763_rast <- rasterize(vect(mystudyarea3763), mystudyarea3763_r, touches=TRUE, field = 1)

# Resample IDW to match the template
idw_r <- resample(idw_r, mystudyarea3763_rast)

# Mask: keeps values only inside the study area, 
# but ensures the full extent is covered
idw_r <- mask(idw_r, mystudyarea3763_rast)
plot(log(idw_r+1))

(r5 <- ggplot(mystudyarea3763) +
    geom_sf(fill=NA, col="grey") +
    theme_void() +
    geom_spatraster(data = log(idw_r+1)) +
    geom_sf(fill=NA, col="black", lwd=2) +
    geom_sf(fill=NA, col="white") +
    scale_fill_viridis_c(na.value = NA))

(r5b <- r5 + 
  geom_sf(data = occ_sf3763, colour = "white", fill = NA))




# Kernel density
library(spatialEco)

crs(myr10k_) <- crs(occ_sf3763)

pt.kde <- sf.kde(x = occ_sf3763, 
                 bw = 20000, 
                 standardize = TRUE, 
                 res=1000)

pt.kde <- mask(pt.kde, mystudyarea3763)


plot(pt.kde)

(r6 <- ggplot(mystudyarea3763) +
    geom_sf(fill=NA, col="grey") +
    theme_void() +
    geom_spatraster(data = log(pt.kde+1)) +
    geom_sf(fill=NA, col="black", lwd=2) +
    geom_sf(fill=NA, col="white") +
    scale_fill_viridis_c(na.value = NA))

(rbb <- r6 + 
    geom_sf(data = occ_sf3763, colour = "white", fill = NA))